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Efficient  numerical  simulations  of  microstructure  development  in  magnet orheolog- 
ical  (MR)  fluids  are  conducted.  The  simulations,  which  are  based  upon  a  fast 
multipole  algorithm ,  treat  the  magnetic  inclusions  as  two-dimensional  continuum 
magnetic  entities.  The  development  of  microstructure  is  quantified  by  computing 
and  recording  the  time  evolution  of  the  effective  permeability  of  the  composite 
fluid.  Such  a  principle  has  been  previously  exploited  for  the  experimental  mea¬ 
surements  of  microstructure  development  [Jolly,  Bender  and  Mathers,  ERMR’97, 
Yonezawa,  Japan  1997].  As  was  observed  experimentally,  numerical  simulations 
reveal  the  evolution  of  microstructure  to  be  multimodal  in  nature.  Unlike  the  ex¬ 
periments,  the  numerical  simulations  afford  us  the  ability  to  observe  the  physical 
mechanisms  associated  with  various  modes. 


1  Introduction 


The  nature  of  microstructure  formation  in  controllable  fluids  has  been  a  topic  of 
recent  considerable  interest.  It  is  generally  believed  that  the  field  responsive  rheo¬ 
logical  effect  is  sensitive  to  the  nature  of  the  microstructure1.  Further,  it  has  been 
postulated  that  the  time  associated  with  structure  formation  is  an  important  con¬ 
stituent  of  the  overall  time  required  for  the  development  of  stress  within  controllable 
fluids2.  There  has  been  considerable  work  to  understand  the  nature  and  temporal 
behavior  of  structure  formation  in  elect  rorhcological  (ER)  fluids  through  analyti¬ 
cal  models3,4’5  and  numerical  simulations6,7  .  Fewer  analogous  magnetorheological 
(MR)  fluid  studies  are  found  in  the  literature8,9. 

In  a  previous  paper10,  we  have  measured  the  time  scales  associated  with  mi¬ 
crostructure  development  in  MR  fluids  by  examining  the  evolution  of  fluid  perme¬ 
ability  (via.  polarization  measurements)  in  response  to  a  step  change  in  applied 
magnetic  field.  This  experimental  technique  is  somewhat  analogous  to  the  use  of 
dielectric  measurements11,12,13  to  infer  microstructure  development  in  ER  fluids  in 
that  the  dielectric  response  of  an  ER  fluid  is  studied  in  response  to  an  electrical 
stimulus.  These  studies  of  MR  fluids  indicated  that  the  microstructure  response 
was  well-fitted  with  a  bi-exponential  function  where  the  two  time  scales  differed 
by  about  a  factor  of  five.  The  first  experimental  time  constant  was  found  to  be 
somewhat  consistent  with  predicted  flocculation  times  based  on  the  computed  time 
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of  collision  of  two  dipoles  in  a  viscous  media,3.  Inconsistency  was  found  in  the  rela¬ 
tionship  to  volume  fraction  <j>  -  the  simple  theory  predicts  t  ~  <j>~n  where  n  =  5/3, 
whereas  n  6  [2/3,  4/3]  was  observed  experimentally. 

The  notion  that  microstructure  formation  occurs  at  multiple  time  scales  has 
been  observed  in  the  work  of  others.  Halsey  and  Toor4  used  some  elegant  theory 
to  argue  that  structure  formation  in  ER  fluids  is  a  two  step  process:  an  initial 
aggregation  to  a  chain-like  structure  resulting  from  electrostatic  interaction  and  a 
phase  separation  to  a  thermodynamic  ground  state.  The  latter  process  is  a  result 
of  thermally  induced  fluctuations  within  the  chains  and  occurs  on  thermal  diffusion 
time  scales.  Through  numerical  simulations  that  neglected  thermal  forces,  Mohebi 
and  colleagues8  showed  that  structure  formation  in  MR  fluids  (<j>  =  10%)  occurs 
at  two  distinct  time  scales.  These  time  scales  correspond  to  an  initial  formation  of 
disparate  chains,  and  then  a  migration  of  chains  into  longer  and  thicker  structures. 
The  first  process  was  seen  to  occur  on  the  millisecond  time  scale  and  the  second 
process  was  observed  to  be  one  to  three  orders  of  magnitude  slower  depending 
upon  the  sample  thickness.  Hass',  using  similar  dynamic  simulations  on  ER  fluids, 
demonstrated  two  initial  time  scales  separated  by  an  order  of  magnitude  which 
he  attributes  to,  first,  pair  formation  and,  second,  percolating  column  growth. 
The  experimental  work  of  Jolly  and  colleagues10  exhibited  time  scales  that  are 
comparable  to  those  found  in  the  latter  two  simulation  studies,  but  could  provide 
no  insight  into  the  mechanisms  associated  with  these  time  scales. 

In  this  paper,  a  novel  two-dimensional  simulation  of  MR  fluid  microstructure 
dynamics  is  presented.  This  simulation  treats  the  particle  inclusions  as  two- 
dimensional  magnetic  entities.  The  time  evolution  is  considered  to  be  magneti¬ 
cally  quasi-static  and  magnetostatic  forces  are  derived  from  the  solution  of  (steady) 
Maxwell’s  equations,  recomputed  at  each  instant  in  time.  For  this  we  use  a  poten¬ 
tial  theoretic  formulation  where  the  boundary  integral  equations  are  solved  with  a 
fast  multipole  method14.  The  simulations  are  then  used  to  explore  the  multimodal 
nature  of  microstructure  development  in  MR  fluids.  These  results  are  compared  to 
previously  reported  experimental  results10. 

2  Equations  of  Motion 

In  this  section,  we  describe  the  equations  of  motion  for  M  permeable  circular  parti¬ 
cles  in  a  rectangular  container  Q  C  1R2  that  is  filled  with  a  non-magnetic 

viscous  fluid.  The  particles  are  randomly  distributed  initially  and  the  external  mag¬ 
netic  field  Ho  is  applied  in  the  vertical  direction.  The  motion  of  the  kth  particle  is 
governed  by  Newton’s  second  law  of  motion.  Particle  interactions  include  magnetic, 
hydrodynamic  and  repelling  forces  while  Brownian  and  inertial  forces  are  neglected. 
By  performing  dimensional  analysis  and  eliminating  small-scale  terms  (for  details, 
see  Ly  el.al.  "’),  we  obtain 

^-  =  F?«  +  jg>  and  x°  =  xt(0),  (1) 

where  x*,,  t,  F/lag,  F/ep,  are  dimensionless  variables  for  position,  time,  magnetic 
force  and  repelling  force,  respectively. 
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The  spatial  scale  is  R,  where  R  the  the  particle  radius  and  the  temporal  scale  is 


- -,  where  D  is  the  Stokes’  drag  coefficient,  p0  is  the  permeability  for  the  carrier 

/'<:"(! 

oil,  and  H o  is  the  magnitude  of  the  applied  held. 

For  our  simulations,  we  shall  assume  that  both  the  particles  and  the  container 
walls  are  “hareF .  To  approximate  this  regime,  we  shall  follow  the  work  of  Klin- 
genberg  et.al.6  and  propose  that  a  “repelling  force ”  acts  on  the  kt.h  particle  as  it 
approaches  others  or  a  wall  of  the  container.  A  simple  model  for  such  a  force  is 
given,  for  instance,  by 


=  -  J2  exp  (  ^ 


\dki\  -  2 R 


)r hi  ~  exp  (-0 


Krni  -  2  R, 


where  r*;  =  t-z - zr-r,  dki  =  dist(12* ,  12| ) ,  and  0  >  0  is  the  repelling  parameter. 

Ix;  —  x*  I 

The  wall  repelling  force  uses  n*,  an  outward  unit  normal  vector  at  a  point  p  on  the 
boundary  of  the  container  12  where  p  is  nearest  to  x*  on  512  and  |d™all|  =  dist(12*,  12). 

The  magnetic  force  on  12*  can  be  calculated  from  the  local  held  H  with  the  aid 
of  the  Maxwell  stress  tensor  crMax  =  /i0[HHr  —  i|H|2<§],  (<5=unit  tensor)  as 


pmag  _ 
r  k  ~ 


•n*  d.S, 


k  ~  jr2p  '' 

ZTTti0N0R  JaUk 

where  n*  is  the  unit  normal  vector  on  512*  ■  An  accurate  estimate  for  the  magnetic 
force  in  (3)  demands  the  continuous  knowledge  of  the  local  magnetic  held  H,  so  that 
Maxwell’s  equations  must  be  resolved  at  each  instant  in  time.  In  our  simulation, 
we  assume  there  are  no  free  currents  in  the  domain,  so  that  the  local  magnetostatic 
held  H  can  be  written  in  terms  of  a  scalar  potential  <F, 

H  =  -V<F.  (4) 

Moreover,  $  is  the  solution  of  the  Laplace’s  equation 

V  •  (/iV$)  =  0,  (5) 

with  highly  oscillatory  coefficients 


in  the  kth  particle, 
in  the  carrier  oil. 


For  our  simulation,  we  assume  pk  ~  2000/ko-  Along  with  equation  (5)  we  require 
the  continuity  of  the  magnetic  potential  <3?  and  of  the  normal  component  of  B.  That 
is,  for  any  k  =  1,2,...,  M , 


lim 

p  — >■  <912  J: 
p  e  12* 


$(p)  = 


lim 

p  — >■  <912* 
P  e  12* 


<9$ 

Pkjp — (p)  = 

onktP 


lim 

p  — *■  <912* 
p  e  12* 

lim  ; 

p  — >  <912* 
P  £  12* 


<9$ 

POJp - (P), 

onk]V 


where  nk  p  is  an  outward  unit  normal  vector  at  p  6  512*  and  12)r  =  12\12*  is  the 
complement  of  12* . 
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3  Integral  Equation  Formulation  and  Boundary  Element 
Discretization 

Although  the  Coefficients  of  the  Laplace’s  equation  (5)  are  rapidly  changing  in  space, 
they  do  remain  constant  in  each  fl*.  Thus  the  overall  potential  can  be  derived  from 
appropriate  charge  densities  supported  on  the  boundaries  of  the  particles.  These 
densities  satisfy  certain  integral  equations  which  are,  in  principle,  amenable  to 
solution  by  finite  (boundary)  element  approximation.  As  we  discussed  in  detail 
elsewhere15,  the  difficulties  associated  with  the  high  computational  cost  of  classical 
boundary  element  approximation  for  this  kind  of  problem  can,  in  fact,  be  overcome 
through  the  implementation  of  the  Fast  Multipole  Method14. 

To  derive  the  integral  equations,  let  us  denote  by  the  boundary  of  the 
domain  Q.  We  also  denote  by  fig  the  fluid  region  and  impose  the 

following  Neumann  boundary  condition  on  5f2o 

<9$ 

-7^\dn0=a-  (9) 

To  guarantee  the  solvability  and  uniqueness  of  the  solution  of  equations  (5)-(9),  we 
make  the  following  requirements  on  g  and  <f>;  namely, 


g  ds  =  0 


$  dx  =  0. 


A  potential  $  satisfying  (5)-(9)  can  be  represented  by  single-layer  potentials16  in 
the  form 


M 

®(P) 

•  n  Jdt 


G{p,q)tj{q)ds(q),  pen. 


Here,  G(p,  q)  =  4^  log  \p  —  q\  is  the  fundamental  solution  of  the  Laplace’s  equation 
in  IR2 .  The  functions  fj’ s  on  are  appropriate  (unknown)  surface  densities. 

Note  that  the  potential  $  in  (10)  automatically  satisfies  A<3?  =  0  on  for  k  = 
0, 1, .  . .  ,M,  and  the  continuity  condition  (7)  at  the  interfaces.  In  addition,  using 
the  jump  relations  of  potential  theory16,  we  obtain  from  equations  (8)  and  (9)  the 
following  system  of  Fredholm  equations  of  the  second  kind, 

M  n  r\ 

£o{p)  -  2  V  /  7p-G(p,  q)£j (q)ds(q)  -  -2 g(jp),  (11) 

i=0daoJ  vnP 

M  n  C\ 

tk{p)  ~  2AS,  ^  G{p,q)£j(q)ds(q )  =  0,  (12) 

j=0JdV.3  0nP 

where  A*  =  — - —  and  equations  (11 )-( 12)  hold  for  p  £  and 

/'.<•  •  /'g  "  "" 

respectively.  Equations  ( 1 1 )-( 12)  are  then  discretized  and  solved  with  the  help 

of  the  fast  multipole  method14.  A  similar  approach  with  details,  implemented  for 
MR  applications,  can  be  found  in  our  earlier  work15,  where  the  Dirichlet  boundary 
conditions  were  assumed  at  the  exterior  boundary  and  the  resulting  formulations 
involved  both  the  single-  and  double-layer  integral  potentials.  We  also  remark  here 
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that  because  our  the  potential  formulation  $  in  (10)  consists  of  only  the  single 
formulation,  numerical  implementation  is  more  efficient  and  integral  singularities 
are  removable.  Numerical  results  are  discussed  in  the  next  section. 

4  Numerical  Results  and  Discussion 


5%  Volume  Fraction 


25%  Volume  Fraction 


large  time 


Figure  1.  Microstructure  formation  for  5%  and  25%  volume-fraction  samples  in  four  stages 
(/  =initial,  first  time  scale  (ti),  second  time  scale  (T2),  and  large  time). 


Simulations  were  conducted  on  six  different  volume  fractions  ranging  from  5%  to 
30%.  The  external  field  is  applied  in  the  vertical  direction.  Simulation  results  for 
two  volume  fractions  responding  to  a  step  field  input  are  shown  in  Figure  1.  Dis¬ 
crete  times  within  the  simulation  are  shown  from  the  initial  configuration  to  the 
near  steady  state  microstructure.  An  issue  in  such  dynamic  simulations  involves  the 
means  by  which  the  state  of  the  microstructure  is  quantified.  Researchers  have  used 
quantities  related  to  nearest  neighbor  inter-particle  distances  and  spatial  correla¬ 
tion  functions  to  monitor  the  state  of  a  microstructure  17’5.  Others  have  used  mean 
square  particle  displacement6  '.  In  a  manner  analogous  to  a  previously  reported 
experimental  technique,  we  will  use  the  composite  magnetic  permeability  to  quan- 
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Figure  2.  Simulated  effective  permeability  fi(t)  as  a  function  of  non-dimensional  time  t  for  six 
different  volume  fractions. 


tify  the  evolution  of  microstructure.  This  method  exploits  the  fact  that  a  percolated 
microstructure  is  significantly  more  permeable  than  a  randomly  mixed  microstruc¬ 
ture.  Some  theoretical  limits  of  this  behavior  are  discussed  in  Simon  et.al.18.  It 
was  reported  by  Ly  et.al.15,  that  the  evolution  of  permeability  in  microstructure 
formation  is  closely  related  to  the  mean  particle  velocity.  During  the  simulation 
process,  we  compute  the  effective  permeability  which  represents  the  overall 

response  for  the  microstructure  formation.  The  definition  and  the  formulation  for 
the  effective  permeability,  which  is  based  on  the  theory  of  homogenization19,  are 
derived  in  Ly  et.al.20.  Figure  2  shows  the  permeability  response  of  the  six  samples 
to  a  step  input  in  magnetic  field.  It  should  be  pointed  out  that  the  simulations  are 
displayed  in  dimensionless  times  and  dimensional  times  can  be  obtained  by  multi¬ 
plying  the  dimensionless  times  with  the  time  scale  described  earlier  in  Section  2.  As 
expected,  both  the  initial  and  the  steady  state  permeabilities  are  linearly  related 
to  the  volume  fraction,  <j). 

As  was  found  experimentally,  the  permeability  response  is  well-fitted  with  a 
bi-exponential  function.  These  fits  are  also  applied  to  the  permeabilities  shown  in 
Figure  2  and  the  resultant  time  constants  and  coefficients  are  presented  in  Table  1. 
Microstructure  frames  corresponding  to  the  two  time  constants  (ti  and  t 2)  are 
shown  in  Fig.  1.  From  these  frames,  it  is  evident  that  the  first  mode  corresponds 
to  particle  pair  and  short  chain  formation  and  the  second  mode  corresponds  to 
coalescence  of  the  short  chains  into  longer  percolating  chain  structures.  Similar 
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observations  from  simulations  of  ER  and  MR  fluids  have  been  made  by  others 
including  Hass'  and  Mohebi  et.al.8.  The  two  time  constants  are  separated  by  a 
factor  of  5  to  10.  An  additional  observation  can  be  made  in  noting  that  magnetic 
energy  is  proportional  to  permeability  for  linear  materials.  In  particular,  it  can 
be  seen  that  about  65  —  70%  of  the  microstructures  magnetic  potential  energy 
is  stored  upon  initial  application  of  the  field.  The  balance  of  the  stored  magnetic 
energy  occurs  during  the  microstructure  formation.  For  high  volume  fraction  fluids, 
over  90%  of  the  energy  is  stored  within  the  first  time  constant  dropping  to  below 
70%  as  volume  fraction  decreases. 


fi(t)  : 

Bi- Exponential  Fit 

=  A  *  (1  -  e1/'  )  +  B  *  (1  —  e*/r2)  +  /,(0) 

n 

72 

A 

B 

MO) 

5% 

16.2 

88.8 

0.170 

0.368 

1.46 

10% 

8.81 

43.8 

0.467 

0.667 

1.97 

15% 

5.87 

46.1 

1.090 

0.517 

2.54 

20% 

3.93 

44.2 

1.550 

0.418 

3.22 

25% 

3.35 

38.6 

1.960 

0.281 

3.95 

30% 

3.12 

29.9 

2.030 

0.359 

4.83 

Table  1.  Effective  permeability  and  hi  -exponent  ial  fit. 


Figure  3  depicts  the  first  time  constant  as  a  function  of  volume  fraction  <f>. 
Simulation  results  as  well  as  experimental  results  from  Jolly  et.al.10  are  presented 
where  the  simulation  results  have  been  dimensionalized  to  correspond  with  the  ex¬ 
perimental  conditions.  As  shown  in  several  references'5  2  0 ,  simple  theory  suggests 

that  the  relationship  between  microstructure  formation  time  and  volume  fraction 
should  follow  a  power  law  behavior  with  a  power  index  of  n  =  —5/3.  In  particular, 
it  has  been  argued  that  the  time  for  pair  formation  is  proportional  to 

“~f ■I#1'* <13> 

where  Jv  is  the  particle  polarization.  Our  simulation  results  have  indicated  that 
the  first  time  constant  is  inversely  proportional  to  volume  fraction  ( n  =  —1)  and 
experimental  results  have  yielded  power  indices  between  —4/3  and  —2/3.  The 
disparity  of  this  behavior  with  the  simple  theory  is  not  surprising  since  the  theory 
is  based  on  two  isolated  particles  (magnetically  modeled  as  dipoles)  in  a  viscous 
medium.  We  further  note  that  our  simulations  do  not  account  for  the  magnetic 
nonlinearity  of  the  particle  material,  which  may  contribute  to  the  range  of  power 
law  indices  in  the  experimental  data.  If  it  is  assumed  that  Jp  =  B,;/</>,  where  B,;  is 
the  intrinsic  induction  of  the  composite,  then  B,;  can  be  substituted  into  Eq.  (13) 
and  the  aforementioned  power  law  index  becomes  n  +  2  (hence,  ii  i  2  -  -  1/3  for 
theory,  n  +2  =  1  for  simulation  and  2/3  <  (n  +  2)  <  4/3  for  experiments).  Figure  4 
contains  previously  reported10  experimental  time  constants  and  scaled  simulation 
time  constants  as  a  function  of  r)4>/Bf.  It  can  be  seen  that  the  time  constants 
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rise  linearly  with  this  function.  The  numerical  results  are  shifted  with  respect  to 
experimental  results.  This  may  be  the  result  of  structured  error  in  the  experimental 
measurement  of  S,; . 


Conclusions 

A  numerical  algorithm  has  been  developed  to  accurately  study  both  the  microstruc- 
ture  and  the  permeability  with  the  aid  of  the  fast  multipole  method.  The  results 
have  been  quantified  with  the  two  time  scales  for  structure  formation  (iq;  particle 
pairing  and  short-chain  formation;  72:  coalescence  of  short  chains  into  longer  per¬ 
colating  chain  structures).  Our  investigation  reveals  that  the  first  time  scale,  Tf,  is 
inversely  proportional  to  the  particle  volume  fraction  -  a  result  that  is  in  reasonable 
agreement  with  previous  experimental  results10.  It  is  anticipated  that  the  use  of 
permeability  (or  polarization)  measurements  to  monitor  the  state  of  microstruc- 
ture  in  MR  fluids  will  be  particularly  useful  in  shear.  Other  means  of  quantifying 
microstructure,  especially  those  involving  monitoring  mean  particle  motion,  may 
fall  short  in  shear  environment.  Permeability  measurements  from  both  particle 
simulation  and  experiments  should  be  capable  of  resolving  microstructure  forma¬ 
tion  and  the  shear-induced  microstructure  degradation.  Correlation  between  such 
measurements  and  the  MR  response  will  be  a  topic  of  future  research. 


Figure  3.  The  first  time  constant  as  a  function  of  volume  fraction:  (triangles) -experimental results 
for  H  —  17.5  kA/m  and  rj  =  0.13  Pa-s;  (diamonds) -experimental  results  for  H  =  7  kA/m  and 
7]  =  0.13  Pa-s.  The  small  circles  are  simulation  results  that  have  been  dimensionalized  to  match 
the  corresponding  experimental  results.  Inverse  linear  fits  are  also  shown.  Experimental  results 
are  from  Jolly  et.al10. 
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Figure  4.  Microstructural  response  times  as  a  function  of  their  assumed  theoretical  dependency. 
Triangles  and  diamonds  are  the  first  and  second  dimensionalized  time  constants  from  simulation. 
Circles  and  dots  are  the  first  and  second  experimentally  measured  time  constants10. 
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